Supersonic crack propagation in a class of lattice models of Mode III brittle fracture 
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We study a lattice model for mode III crack propagation in brittle materials in a stripe geometry 
at constant applied stretching. Stiffening of the material at large deformation produces supersonic 
crack propagation. For large stretching the propagation is guided by well developed soliton waves. 
For low stretching, the crack-tip velocity has a universal dependence on stretching that can be 
obtained using a simple geometrical argument. 
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The determination of crack-tip velocities in brittle 
fracture is one of the unsolved problems in fracture 
mechanics Q, Q- Spontaneous crack propagation typi- 
cally occurs when some energy balance condition is sat- 
isfied. This is usually known as the "Griffith's criterion" 
and it states roughly speaking, that the available elas- 
tic energy has to be larger than the energy necessary to 
create new surface in the material upon crack advance. 
We will define the excess elastic energy (EEE) as the 
difference between the available elastic energy and the 
(static) crack-surface energy. When the EEE is positive 
a dynamic regime is achieved. The nature of this regime 
depends on the value of the EEE [J H 0. If the EEE 
is low, a stationary regime is usually observed, in which 
the crack velocity is a well defined (material dependent) 
function of the EEE, and if the geometric configuration 
is symmetric with respect to the crack, crack propaga- 
tion occurs along a straight line. For larger values of the 
EEE, crack path oscillations, crack branching and other 
phenomena typically occur. 

We will focus here in the low EEE regime. Even in 
this relatively simple case, the prediction of crack-tip ve- 
locity as a function of material parameters, sample ge- 
ometry and strain conditions is a challenging problem. 
Linear elastic fracture mechanics (LEFM)yJ,|2j suggests 
that the crack-tip velocity V should be a univocus func- 
tion of the EEE and material characteristics. Among 
these, microscopic details of the process zone must be in- 
cluded in general. Another most important prediction of 
LEFM is the impossibility that the value of V exceeds the 
Rayleigh velocity of the material, which is typically some- 
what lower than the shear sound velocity. There is some 
experimental evidence that challenges thispredictionQ. 
In addition, it has been recently shown|a| that elastic 
stiffening at large deformation can unambiguously pro- 
duce cracks running faster than the shear sound velocity 
in simulations of mode I (opening) propagation. 

Here we analyze the dependence of crack tip velocity 
upon model parameters in a class of piece-wise harmonic 
lattice model for mode III fracture. We find supersonic 
propagation and study its characteristics. We find a re- 
markable regime in which the crack-tip velocity can be 
accurately predicted on the basis of a geometric argu- 
ment, and is independent of most of the model parame- 



ters. 

We model a stripe of material, laterally stretched in the 
out-of-plane direction, with a crack propagating along 
the mid line. This geometry is particularly adequate to 
study stationary crack propagation as the system looks 
identical upon extension of the crack. This kind of mod- 
els have been studied before, for instance in Refs. 0,0- 
The model consists of a discrete square lattice. At each 
node a mass is placed, and the links to nearest neigh- 
bor sites are modeled by piece-wise linear springs. The 
full lattice forms a rectangular system, to which out-of- 
plane fixed-displacement lateral boundary conditions are 
applied. We take the horizontal x axis along the stripe. 
This is also the crack propagation direction. The out- 
of-plane displacement of the masses (the only one which 
is non-zero) is noted by u. The equation of motion of a 
mass at the discrete position i,j in the lattice is simply 
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where [i',f] are the neighbor sites to i,j and the func- 
tion F specifies the force law of the corresponding spring: 
F(x) — x for x < u n i, and F(x) = jx for x > u n i for in- 
tact springs, and F{x) = for broken springs. Note that 
the low deformation spring constant and mass of the par- 
ticles have been set to 1 (this sets the wave velocity to 1, 
also), and that u n i is the critical stretching at which the 
spring constant changes from 1 to 7 (7 > 1). 

Springs are assumed to break irreversibly when the 
stretching is larger than a fixed threshold value Ubk- To 
avoid crack branching, which may occur under certain 
conditions, we allow to break only the vertical springs 
in the center of the stripe (for low stretching this pro- 
cedure is not necessary, see below). For symmetry rea- 
sons we simulate only the upper half of the stripe, and 
N y will denote the number of rows in this half sys- 
tem. The boundary condition at the bottom border 
is Uii — —Uifi, whereas at the upper border we have 
Ui,N y +i = (N y + 1/2)6, where 6 is the nominal stretching 
applied to each vertical spring. We always keep 6 < u n i, 
not to force the system in the non-linear regime with the 
uniform applied strain. Along the x direction, we take 
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the boundary conditions as follows: ujv T .j = {j — 1/2)5 
and uoj = (N y + 1/2)5. These conditions correspond to 
those that should be valid at x — ► ±oo. By enforcing 
them at a finite cc-position some error is introduced, but 
we check that this error is small if N x is large enough. 
The crack tip is initially placed at an x position close 
to the right border of the system, by seeding the simu- 
lations with an appropriate initial condition. Along the 
simulation the crack advances in the positive x direc- 
tion. Each time we detect that a new spring ahead of 
the crack has broken, we shift all system coordinates and 
velocities one unit to the left. In this way the crack tip 
is always maintained at the same position with respect 
to the lattice, and the instantaneous crack-tip velocity is 
determined as the inverse of the difference in consecutive 
times of spring breaking. Reported values correspond to 
averaging over long simulation times, after dependences 
on initial conditions have died out. 

In the present stripe geometry, the Griffith's crite- 
rion is stated as follows: crack propagation cannot oc- 
cur if stretching is lower than some critical value 8q- 
Equating the elastic energy far ahead of the crack-tip, to 
that far behind it, we find 8g — v2TV (2-^y + 1); where 
T = ul k /2 + (7 - l)(ubk - u n i) 2 /2 is the crack energy 
per bond. We see that S G -> for N y -> 00. The EEE 
per horizontal bond (referred to as e) is easily calculated 
in terms of 8g as e = T [(S/Sq) 2 ~ l] • As it should, 
this quantity vanishes at 5 = Sq. Note that contrary to 
the usual assumption in LEFM, in our model the excess 
energy e is not dissipated at the crack tip, but it is trans- 
formed into kinetic energy, that eventually escapes from 
the left border of the system. 

To test the classical prediction in the normal (sub- 
sonic) case, in Fig. ^ we show curves for the velocity 
as a function of stretching for 7 = 1 and different values 
of N y . The curves converge to a well defined behavior 
for large N y if the plot is done as a function of e. The 
result agrees with the predictions of LEFM. The velocity 
approaches the wave speed for large e, whereas at low e 
is reduced due to lattice trapping effects. In accordance 
also with previous simulations in this kind of models)?]], 
the instantaneous velocity (Fig. [I] inset) is found to tend 
to a constant for large stretching, but it is a fluctuating 
function of time in the low stretching regime. 

Results change if 7 > 1. In this case different velocity 
curves for different stripe widths do not tend to collapse 
in the large N y limit if plotted as a function of e. Instead, 
they do if plotted as a function of the unperturbed value 
5 of the strain ahead of the crack. This is an important 
result: since the propagation is supersonic, the crack can 
respond only to the local state of the medium, and the 
velocity must be independent of the system width if this 
is large enough. In Fig. 0we plot the velocity dependence 
on S for different system parameters. We clearly see the 
supersonic propagation, which occurs for any 5 > if N y 
is large enough. 
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FIG. 1: (Color online) Average velocity of the crack in the 
harmonic case as a function of the excess elastic energy per 
horizontal bond e for different stripe widths N y . To check for 
independence on stripe length, two different values of N x are 
reported. In the inset, values of the instantaneous velocity as 
a function of time for a system of N x = 320, N y = 160 are 
indicated in two cases. 



1 .6 F ^ 

. 1.05 
1.00 



1.4 













^5=0.3 












,j>=0.2 




0.00 0.25 0.50 0.75 1.00 

8/1/, 

FIG. 2: (Color online) Velocity of the crack in a system with 
N x = 120, N y — 240 for different model parameters as a func- 
tion of stretching. The thick continuous line is an analytical 
fitting valid for low 8, large 7, and N y — » 00. The inset shows 
that the fitting improves by increasing system size (N y —* 00 
corresponds to V*(6g) —> 1, the values predicted by the fitting 
are indicated by the arrows) . The points reported in the inset 
correspond to all combinations of the parameters: 7 = 2,4, 
N v = 150,240,360, and u bk = 1.5,2,2.5,3, plus those in the 



main plot. 



To analyze these results, it is convenient to solve first 
the following formal problem. Consider a perfectly har- 
monic, infinite and continuous system, stretched some 
fixed amount S per unit length in the y direction. Sup- 
pose that by some external mean we force a crack to prop- 
agate along the x direction at some velocity V, which is 
arbitrary, with the only condition V > 1. The displace- 
ment field u(x, y) is found to be simply: 



u(x,y) = yS 



outside the Mach's cone 
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u(x,y) = ± = inside the Mach's cone, for y ^(0) 

V V 2 — 1 

Here, the crack tip is located at x = 0, y = 0, and the 
Mach's cone half-angle 9 is given by sin^ = 1/V. Note 
that inside the Mach's cone this solution has a constant 
value of \du/dx\, given by \du/dx\ = S/^/V 2 — 1. In 
the case in which the same problem is solved (numeri- 
cally) on a discrete square lattice (inducing the propa- 
gation by breaking bonds at a fixed rate V), an oscil- 
lation close to the border of the Mach's cone is found, 
and the maximum stretching |Au| m ax of a spring in the 
system is enhanced by a constant factor. It is found to 
be | Am | ma x ~ 1.26 S/VV 2 - 1. Note that |Azi| max di- 
verges when V —> 1 no matter how small S is. 

In Fig. |31 we see the actual results from the simula- 
tions in the non-linear model. At low 5 (upper panel) 
the Ut,j function looks roughly similar to that given by 
Eq. J2| (except for the reflection effects at the lateral 
borders of the system). We find a number of spatial 
positions at which the non-linear threshold is exceeded. 
This is quite reasonable, otherwise we should not expect 
any supersonic propagation. In addition, it is found that 
the instantaneous crack tip velocity is a fluctuating func- 
tion of time, typically around ± 10 % of the mean value. 
This fluctuation is acoompanied by a change in the lo- 
cation of the points at which the non-linear threshold is 
exceeded. In spite of this, the mean velocity can be deter- 
mined to a high precision by averaging over a long sim- 
ulation time. Upon increasing S, the non-linear regions 
tend to arrange in the form of solitons, that eventually 
(for S/u n i between 0.7 and 0.75 in Fig. |3J) go outside 
the Mach's cone. There are five solitons at S/u n i = 0.75 
in Fig. |2| and this number reduces as 5 increases. For 
the largest S values studied, typically a single soliton is 
observed, which drives the supersonic propagation. We 
have carefully verified that the solitons in front of the 
crack tip are simply the non-linear solitary waves corre- 
sponding to the present non linear elastic system, i.e., 
they are the analogous of the solitons appearing for in- 
stance in the well known Toda lattice 8] . In the regime 
of solitons-driven propagation the instantaneous veloc- 
ity is constant. The transition between the low-<5 and 
the soliton driven regimes is found to be abrupt, with 
some hysteresis upon increasing/decreasing 5, and with 
a small but clearly observable jump in the velocity (not 
appreciable in the scale of Fig. 0) . 

The low 5 region of the plot in Fig. shows a re- 
markable independence on model parameters. Actually, 
this part of the curve can be understood and fitted on 
the basis of the following heuristic argument. Consider 
the previously found solution for the harmonic system, 
Eq. J5J). We look for the possibility that this solution 
is self-maintained in the presence of non-linearities (i.e., 
when 7 / 1), and we analyze in particular the case in 
which 7 is very large. It is then found that the velocity 
V* obtained by requiring u n i = |Au| m ax plays a spe- 



cial role. This velocity is (taking into account the factor 
introduced by the discreteness of the system) 

V*{5) ~ v/l + (1.26 S/u nl ) 2 . (3) 

In fact, if V > V*, the solution does not explore the non- 
linear part of the potential anywhere, and then it can- 
not be self- maintained. On the other hand, if V < V*, 
there are regions in the system in which the non-linear 
threshold is exceeded by a finite amount, but this cannot 
be acceptable if 7 is sufficiently large. Then we expect 
that the velocity at which the crack propagation stabi- 
lizes is precisely V * . This prediction is plotted on top of 
the numerical results in Fig. [21 The only parameter of 
the model on which V* depends upon is the non-lineal 
threshold u n i , beyond that the solution is independent of 
the precise values of 7 (assumed large) and Ubk ■ The fit- 
ting improves for larger systems as shown in the inset to 
Fig. |3 where different points for all combinations of the 
parameters indicated are plotted as a function of V*(5g)- 
The continuous lines in Fig. [21 (inset) correspond to the 
finite size ansatz V(5) = V*(S) - V*(5 G ) + 1. The nu- 
merical data follow accurately this trend, and numerical 
extrapolation for N y — > 00 allows to claim that the fit- 
ting J2Jl is better than 1 % for S/u n i < 0.3 for all the 
parameters studied. This accuracy is remarkably good 
because, as we already said, the instantaneous velocity 
in the low S regime fluctuates around the mean value in 
about 10 %. 

An attempt to predict by the same kind of geomet- 
ric argument the propagation velocity of intersonic mode 
I or mode II cracks gives negative results: A theoreti- 
cal investigation in a purely harmonic system shows that 
(contrary to the mode III solution J2J) the stress fields di- 
verge at the crack tip and are dependent on the boundary 
conditions at infinity (see !2|, pp. 348-355). Then there is 
not a well defined maximum stress to be used in an argu- 
ment similar to that presented for mode III. In addition, 
the consideration of truly supersonic mode I or mode II 
cracks for totally harmonic springs up to breaking shows 
that (in opposition to the mode III case) for any fixed 
velocity all stresses fall below any chosen threshold value 
if S is taken small enough, and this implies that super- 
sonic propagation cannot occur for arbitrarily low values 
of stretching. These considerations show that the present 
mode III problem stands as a very particular but impor- 
tant case in which the crack tip velocity can be predicted 
on the basis of general arguments. 

The simulations that produced Fig. El were done by 
allowing to break only the vertical springs ahead of the 
crack. However, we can check a posteriori whether other 
springs overpass the breaking threshold Ubk or not. It is 
found in general that there is a separating value So below 
which the previous prescription is not necessary. If all 
springs are allowed to break, there are no changes in the 
results for 5 < 80, whereas other phenomena (typically 
crack branching) are observed if 6 > 6q. We found Sq 
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FIG. 3: (Color online) Four views at different S values of the 
supersonic cracks, for 7 = 4 and Ubk/u-ni — 2 in a system 
of size 120 x 240. The mj function is plotted (for clarity, 
we plot the mesh using only 1 every 6 x 12 actual points 
in the left halves of the plots), and the surface is shadowed 
according to the values of (m+ij — Uij)/u n i. The points of 
the system at which the non-linear threshold u„i is exceeded 
(i.e, when (ui+ij — Ui t j)/u n i > 1) are highlighted in black on 
the right halves of the plots. The configuration of the upper 
two panels are not totally stable. Instantaneous velocity as 
a function of time fluctuates as shown for 8/u„i = 0.3, and 
the spatial location of points where the nondinear threshold 
is exceeded changes with time. 



the non-linear threshold for the vertical springs is moved 
to infinity, the velocity curves obtained are only slightly 
modified, the low 8 fitting Jjyl remains good, and in par- 
ticular propagation remains supersonic. On the other 
hand, if the non-linear threshold of the horizontal springs 
is moved to infinity, supersonic propagation completely 
disappears. Then we arrive at the seemingly paradoxical 
result that the non-linearities of the springs that actu- 
ally break are not crucial in determining the propagation 
velocity. This consideration is important for theoretical 
arguments since it tells that the present results cannot 
be explained with an analysis assuming a Barenblatt- 
type process zone (see 0, ch. 3), since in this case only 
non-linearities in the vertical links are considered. 

In conclusion, we have studied the supersonic propaga- 
tion of cracks in a lattice model of mode III fracture, in 
the context of elastic stiffening at large deformation[{|. 
Crack velocity is found to depend on the local strain 
ahead of the crack. For large stretching, the crack prop- 
agation is driven by solitons formed in the non-linear lat- 
tice. In the low stretching regime well developed driving 
solitons are absent. In this last case the crack velocity 
can be accurately predicted on the basis of a geometri- 
cal argument, and it is found to have a general explicit 
expression. This stands as one of very few predictions of 
crack propagation velocities in models of brittle fracture. 

The authors acknowledge financial support from CON- 
ICET (Argentina). 



is roughly 0.2 ~ 0.3 for the parameters we have simu- 
lated. As supersonic propagation occurs in our model for 
any S (in the N y — > 00 limit) we have here an example 
of supersonic crack propagation in which the crack tip 
is stable without the ad hoc introduction of breakable 
springs located only on a previously defined crack path. 

All simulations presented have been done in the ab- 
sence of any dissipative term. Therefore we can clearly 
ascribe the supersonic crack propagation to the non- 
linearities of the potential, in contrast with other models 
in which supersonic propagation is observed to be in- 
duced by the existence of some kind of dissipative terms 
. If in our model a dissipative Kelvin termy} of typical 
strength a is included (i.e., a generic term at the r.h.s. 
of Eq. of the form ~ ad(W 2 u)/dt), we have verified 
that the behavior of the velocity is smooth when a — > 0, 
re-obtaining the results with no dissipation in this limit. 

It is absolutely clear that the stiffening of the potential 
is the responsible for the supersonic crack propagation 
in this model. However, we have observed that the ef- 
fect of the stiffening of horizontal and vertical springs on 
the propagation is very different, and in a certain sense, 
counter intuitive. Non-linearities of vertical springs do 
not have a qualitative effect on the results presented: If 
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